Displacement Discontinuity Research Articles

Splitting proximate algorithm for deformable image registration based on functions of bounded deformation.

Deformable image registration is a crucial task in the field of medical image analysis. Functions of bounded deformation (BD) have been proved effective for modeling the displacement fields between medical images since they can capture the discontinuity of displacement fields along edges of organs andtissues. However, we find that at the same time, BD functions-based models tend to obtain discontinuous displacement fields inside the regions of organs and tissues due to image noises in some cases and the presented gradient descent algorithm is time-consuming. To alleviate these problems, we propose a faster algorithm named SPA: splitting proximatealgorithm. In the framework of variable-splitting scheme, we incorporate a proximal term in the deformable registration energy based on functions ofBD. The proposed algorithm can efficiently solve the original model and obtain displacement fields, which look more natural and plausible. Numerical experiments show the effectiveness and stability of the proposed algorithm. The proposed SPA is able to drastically register the images with a plausible deformation field and not sensitive tonoise.

Modification of Segment Structure Calculation Theory and Development and Application of Integrated Software for a Shield Tunnel

From the aspect of calculation theory, the beam–spring model method and modified routine method of shield tunnel segment structure calculation were improved, and an efficient integrated software system for segment structure calculation of shield tunnel was developed. The beam–spring method is generally calculated according to the assumption of continuous displacement between beams and joints, and the existing modified routine method assumes that the lateral pressure gradient is constant generally, which does not consider the variation in lateral pressure gradient caused by the difference in the lateral pressure coefficient of soil layers or the water level height, which has a certain deviation from the actual situation. The existing beam–spring method and modified routine method theory were improved, the discontinuous displacement between beams and joints in the beam–spring method was taken into account, and the problem of lateral pressure gradient change in the modified routine method was solved. The calculation software system developed by C# and python programming language was proposed to improve the accuracy and efficiency of segment structure calculation. Based on the actual monitoring data of the internal force of the shield tunnel segment and the adjacent shield tunnel segments under construction in Changsha, China, the segments of the shield tunnel with different cross-section sizes and different hydrogeological conditions are calculated to verify the reliability of the calculation software system. At the same time, combined with the calculation results of the software system and field test data, the stiffness reduction coefficient and equivalent foundation resistance coefficient in the modified routine method were derived to further improve the accuracy of the calculation results, which provided a new idea for the calculation of segment structure of shield tunnel with different diameters under different hydrogeological conditions.

Numerical investigation of the fracture network morphology in multi-cluster hydraulic fracturing of horizontal wells: A DDM-FVM study

Hydraulic fracturing is the most effective method to enhance the recovery of unconventional oil and gas via the formation of complex connected fracture networks. In this study, a pseudo three-dimensional fully fluid-solid coupled model is established to investigate the interaction between hydraulic fractures (HFs) and natural fractures (NFs), and the fracture network propagation process. The displacement discontinuity method (DDM) and finite volume method (FVM) are applied to simulate the rock deformation of the reservoir and fracturing fluid flow into the fracture networks, respectively. The influence of the stress shadow, flow rate distribution, comprehensive fluid loss, and complex extension behavior are considered in this model. A modified analytical crossing criterion is used to predict the interaction between the HF and NF. To improve the accuracy of the pressure calculation, the width of the fracture vertical cross-section is replaced by an equivalent width. The coupled equations are solved by Newton–Raphson iteration method. The numerical results of the proposed model agreed well with published analytical and experimental results. Based on the model, detailed sensitivity analyses are conducted to investigate the interaction between HFs and NFs, and the geometries of the fracture networks. The simulation results show that the fracture network morphology is significantly affected by NFs in naturally fractured reservoirs. The stress interference between HFs can weaken the effect of the horizontal stress difference; thus, NFs can be activated effectively. In addition, the complexity of the fracture network is proportional to the activation rate of the NFs. A more complex fracture network can be formed under the following conditions: a small horizontal stress difference, large NF length, low injection rate and low viscosity. It is noteworthy that a smaller value is not always better for cluster spacing. An optimal cluster spacing exist to obtain good results for reservoir stimulation.

Geomechanical Model for Frictional Contacting and Intersecting Fracture Networks: An Improved 3D Displacement Discontinuity Method

Summary The displacement discontinuity method (DDM) is widely used in large-scale engineering problems, such as hydraulic fracturing stimulation in unconventional reservoirs and enhanced geothermal systems, due to its convenient calculation and high precision. Although the 3D-DDM based on triangular elements can more accurately describe intrinsically geometric characteristics of artificial/natural fractures, due to the disadvantage of constant DDM, it will still produce large errors when calculating frictional contacting, close-spacing, and intersecting fractures, which greatly limits the universality of DDM in complex fracture networks. To the best of the authors’ knowledge, few DDM-based models are capable of tackling partially contacting fracture networks with arbitrary intersecting angles. In this paper, we propose a more efficient 3D-DDM algorithm via integrating the analytical solution, the 20-point Gaussian quadrature formula for standard triangles (GQSTS) integration algorithm, and the adaptive Gaussian-Kronrod integration algorithm. Then, combining the “local mesh refinement” grid, the Mohr-Coulomb correction for negative fracture aperture, and two rough fracture deletion strategies, the optimized 3D-DDM algorithm with broader versatility is established. All optimization measures are validated by the relevant fracture model. In the optimized 3D-DDM algorithm, (1) the “local mesh refinement” technique effectively improves the calculation accuracy of intersecting fractures; (2) the adaptive Gaussian-Kronrod integration algorithm not only improves the integration accuracy of high oscillation functions but is also 50% faster than the traditional Gaussian integration algorithm; (3) the Mohr-Coulomb criterion serves as a rigorous constraint for frictional contacting fractures to avoid the occurrence of negative fracture aperture; and (4) using the rough “deep deletion” or “shallow deletion” strategy to replace the ultrahigh-quality mesh generation can more conveniently maintain the accuracy of the complex fractures model. This optimized 3D-DDM algorithm may serve as a basis for simulating nonplanar 3D fracture swarms’ evolution, allowing for frictional contacting, closespacing, and intersecting fractures.

Displacement discontinuity method for interfacial cracks in one-dimensional hexagonal quasi-crystal coating under thermal-mechanical loading

This study derived Green’s functions of the concentrated interfacial displacement and temperature discontinuities. Green’s functions and the displacement and temperature discontinuity method were used to analyze interfacial cracks in a one-dimensional (1D) hexagonal quasi-crystal (QC) coating under thermal-mechanical loading. The boundary integral-differential equations were established for an interfacial crack, the fundamental solutions were obtained for uniformly distributed temperature, and the displacement discontinuities were applied over a line element. The near-crack-tip oscillatory singularity of stresses was eliminated by introducing the Gaussian distribution function to approximate the delta function. After this treatment, the expressions of the stress and heat flux intensity factors without oscillatory singularity, and the energy release rate (ERR) at the crack tips were obtained. The displacement and temperature discontinuity boundary element method for numerical simulation is proposed, and the influence of relevant factors is briefly discussed.

Numerical Modeling of Fracture Height Propagation in Multilayer Formations Considering the Plastic Zone and Induced Stress.

Hydraulic fracturing and acid fracturing are very effective stimulation technologies and are widely used in unconventional reservoir development. Fracture height, as an essential parameter to describe the geometric size of a fracture, is not only the input parameter of two-dimensional fracturing models but also the output parameter of three-dimensional fracturing models. Accurate prediction of fracture height growth can effectively avoid some risks. For example, petroleum reservoirs produce a large amount of formation water because wrong fracture height prediction leads to the connection between the oil or gas reservoir and the water layer. Although some fracture height prediction models were developed, few models considered the effects of the plastic zone, induced stress, and heterogeneous multilayer formation and its interaction. Therefore, considering the influence of many factors, an improved fracture-equilibrium-height model was developed in this study. The successive over-relaxation iteration method and the displacement discontinuity method were used to solve the model. We investigated the effects of the geological and engineering factors on fracture height growth by using the model, and some important conclusions were obtained. The higher the fracture height, the larger the plastic zone size, and the more obvious its influence on fracture height propagation. High overlying or underlying in situ stress and fracture toughness and low fluid density played a positive role in limiting the growth of the fracture height. Induced stress caused by fracture 1 could not only inhibit the height growth of fracture 2 but also promote its growth. The model established in this paper could be coupled to a fracturing simulator to provide a more reliable fracture height prediction.

A class of strain‐displacement elements in upper bound limit analysis

AbstractThe main restriction in applying the upper bound theorem in limit analysis by means of the finite element method is that once the displacement field has been discretized, the flow rule and the boundary conditions should be satisfied everywhere throughout the discretized structure. This is to secure mathematically that the result will be a rigorous upper bound of the exact limit load multiplier. So far, only the linear and the quadratic displacement elements with straight sides, which result in constant and linear strains, respectively, satisfy this requirement. In this article, it is proven both theoretically and practically that there is a general class of strain‐displacement triangular elements with straight sides which provide rigorous upper bounds. The interpolation scheme is based on the Bernstein polynomials, and there is no upper restriction of the polynomial order. Both continuous and discontinuous displacements can be considered. The efficiency is examined through plane strain and Kirchhoff plate examples. The generalization to 3D (i.e., tetrahedral elements with plane faces) and other structural conditions is straightforward.

Skeleton curve and phantom node method for the Thick Level Set approach to fracture

The Thick Level Set method (TLS) is an approach for non-local damage modeling in which the damage evolution is linked to the movement of a damage front described with the level set method. More recently, a new version of the TLS, designated as the TLSV2, has been proposed as a new concept for coupling continuum damage modeling and discrete cohesive crack modeling for failure analysis in solids. The main objective of this new framework is to profit from both modeling approaches. The continuum part allows for handling crack initiation, branching and merging, whereas the cohesive part brings the capability to handle discrete cracks with large crack opening or sliding without heavily distorted elements, and with the possibility to model stiffness recovery upon contact. In this paper, a generalized framework for the TLSV2 is introduced. Two major issues with the TLSV2 method that have not been dealt with since its inception are addressed in this study, and solutions are proposed. Firstly, the method depends on the location of skeleton curve of the level set field, on which the discontinuity in the displacement field is evaluated. The problem of locating the skeleton curve can be a complicated task, mainly because topological events may emerge as the analysis progresses, such as crack branching. The skeleton curve is determined through a combination of ball-shrinking and graph-based algorithms and then mapped onto the finite element mesh. Secondly, the cohesive forces and displacement discontinuity of the TLSV2 are modeled using the phantom node method. Furthermore, a new approach to compute the averaged values of local quantities is introduced, and model calibration is discussed. The degree of stiffness recovery under compression that is still needed for the continuum part is investigated. Numerical experiments demonstrate the accuracy and ability of the proposed model to handle simulation of failure analysis presenting complex topological crack patterns.

Simulation and optimization of reconstructive surgery procedures on human skin

The complex procedures involved in the reconstructive surgery of human skin to minimize post-operative scarring are here modeled by means of an automated computational tool. A finite strain no-compression membrane model accounting for the tendency to develop wrinkling regions in the skin is presented. The constitutive behavior of the material is then described by a suitable hyperelastic incompressible potential. Transpositions of skin flaps during surgery procedures are here computationally described by a general mapping technique of the internal boundary corresponding to surgery cut. The archetypal reconstructive surgery of a Z-plasty, where a rotational transposition of the resulting triangular flaps is involved, is considered in details, along with multiple Z-plasty and rhombic flap transposition. The results are discussed in terms of optimal deformation parameters, related to stress/strain localization, displacement discontinuities and wrinkling.

Simultaneous numerical simulation of the hydraulic fractures geometry in multi-stage fracturing for horizontal shale gas wells

In this study, the hydraulic fracture geometry of multistage treatments is simulated. A simultaneous numerical simulation strategy is implemented for the three modules that allow for hydraulic fracture geometry calculation. The implemented modules of the simulation strategy are the fracture propagation, the formation stress, and the flow rate distribution, which are coupled to obtain a simultaneous solution. A second-order partial differential equation that describes the fracture propagation module is solved using an implicit Finite Difference Method (FDM), taking advantage of the concept of transmissibility. The flow rate distribution module calculates the inlet flow rate at each fracture by solving a nonlinear equation system. The formation stress module is described by a linear equation system, solved with the Discontinuous Displacement Method (DDM) for induced stress formation calculation. Coupling the modules requires the solution of two implicit equation systems and the evaluation of a fracture length criterion. A case is studied based on a single well (i.e., the X2-HF well) located northeast of the Fuling shale gas field, targeting the Wufeng and Longmaxi formation. A four-stage stimulation treatment with three fracture clusters is simulated. The space discretization for each stage is done for 30 elements with 10-time steps of 666 s. The obtained geometry in the mentioned field has a maximum half-length of 165.6 m and a minimum half-length of 136.8 m. The numerical solution requires 9.55 s to solve a single-stage scenario and 34.2 s to solve a four-stage scenario.

Deformation mechanism of WC single crystals under nanoindentation: Effects of surface defects and orientation on pop-in and hysteresis

Nanoindentation was carried out on pure tungsten carbide (WC) on the basal (0001) and prismatic (1010) planes, using Berkovich and spherical indenters, in both single load and multi-load testing. The work focuses on correlating the load-displacement curves, including elastic to plastic deformation, size effect and hysteresis with the deformation behaviour of WC. With different specimen preparation processes, the elastic to plastic deformation started at different threshold loads: This observation was found to be due to the variation in surface dislocation density. Staircase deformation was observed thought to be caused by dislocation motion and the formation of slip bands; sudden displacement discontinuities in the load-displacement response - associated with dislocation loop nucleation - occurred at, or near the theoretical shear strength. Furthermore, discontinuities in load-displacement curves were also used to confirm that hysteresis loops were a result of plastic deformation, as they were absent when the loading was purely elastic.

A viscoelastic constitutive model for shape memory polymer composites: Micromechanical modeling, numerical implementation and application in 4D printing

In this work, a novel micromechanics-based thermo-viscoelastic constitutive model for shape memory polymer composites (SMPCs) is proposed and applied to four-dimensional (4D) printed SMPCs. The multi-branch constitutive model is used to simulate the time- and temperature-dependent mechanical behavior of the shape memory polymer matrix. According to the elastic–viscoelastic correspondence principle, the equivalent viscoelastic stiffness tensor of the composite is obtained in the micromechanics framework of energy-based effective strain theory and the Mori-Tanaka homogenization scheme, in which a two-parameter interfacial damage model is adopted to consider the displacement discontinuity and traction continuity conditions at the interface between the inclusion and the matrix. The numerical integration scheme for the three-dimensional (3D) viscoelastic constitutive model of SMPCs is presented, and the finite element application is implemented by the user material subroutine UMAT of ABAQUS. After identifying the model parameters with experimental data, the tensile stress-strain curves and stress relaxation phenomena of 4D printed unidirectional SMPCs are successfully described by theoretical simulations. Besides, theoretical simulations also adequately predict the shape memory behavior of 4D printed composites and complex members.

A developed non-smooth coordinate transformation for general bilateral vibro-impact systems.

Discontinuity and non-smoothness of system displacement and velocity caused by mechanical impact make the related research on dynamics of vibro-impact systems very difficult and complex. For the sake of bypassing the problems resulting from impact to some extent, Zhuravlev and Ivanov coordinate transformations were proposed, which can effectively convert the vibro-impact system to one without impact terms. In this paper, a more direct and universal transformation for general bilateral rigid vibro-impact systems is proposed. It is inspired by the main technique of Ivanov transformation, which makes the trajectories remain continuous in an auxiliary phase space. It can be directly applied to common vibro-impact systems, whether the positions of barriers are symmetrical or the restitution coefficients of barriers on both sides are consistent. In particular, this method can also be applied to the unilateral vibro-impact system. Validity of the proposed methodology is examined by means of case studies.

Spatial distribution of mechanical parameters along a fracture interface

Many theoretical and numerical approaches for wave propagation in fractured media use normal and shear fracture specific stiffnesses to represent the complexity of fracture topology during deformation under stress. In this paper, XLiFE++ is used to simulate the interaction between a fracture and elastic wave on a 2D plane, the simulated displacement field around the fracture is converted into traction on the fracture interface, and the jump in displacement on the fracture interface is calculated relative to displacements on both sides of the fracture. Based on the linear slippage model, the specific value between jump displacement and traction is calculated at each fracture point to obtain the fracture stiffnesses and then the stiffness curves along the fracture surface are plotted using the stiffnesses at each fracture point. Using XLiFE++, 144 orthogonally designed parameter combinations are simulated, the simulated displacement field around the fracture with material parameters and fracture stiffnesses is used as a training sample, and a neural network method for estimating the fracture stiffnesses is provided. The verification for the training results indicates that the ANN model can effectively represent the non-linear relationships between Young’s modulus, Poisson’s ratio, excitation frequency of materials, shear stiffness and normal stiffness of the fracture, and displacement field around the fracture and estimate the fracture stiffnesses relative to displacement field around the fracture surface. Therefore, the proposed method is an effective tool to rapidly determine the fracture interface stiffnesses. • Spatial distribution of elastic modulus around fracture •Simulate the interaction between elastic waves and fracture with different stiffness •Calculation of fracture stiffness based on displacement discontinuity theory

Study on fracture initiation pressure of a horizontal well considering the influence of stress shadow in tight sand reservoirs

Multi-stage hydraulic fracturing (HF) is an effective stimulation technique for enhanced hydrocarbon recovery in tight sand reservoirs, and the fracture initiation pressure (FIP) is a key parameter for evaluating the difficulty of HF. Recent studies have shown that the stress shadow caused by pre-existing hydraulic fractures considerably affects fracture initiation. However, few studies have fully analyzed the FIP considering the stress shadow effect. Additionally, factors affecting stress shadow have not been quantitatively evaluated. Therefore, the aim of this study was to establish an efficient model for investigating the impact of stress shadow on the FIP. We developed a comprehensive stress distribution model for the perforated horizontal wellbore and analyzed the fracture initiation mechanism. Subsequently, we applied the displacement discontinuity method (DDM) to simulate the stress shadow generated by hydraulic fractures. Next, we validated the DDM against the analytical solution and analyzed the variation of the FIP along the horizontal wellbore considering the effect of stress shadow. Finally, we evaluated the key parameters affecting the FIP using a well-trained back-propagation artificial neural network. The results showed that stress shadow increased the magnitude of the FIP during the multi-stage HF operation. The fracture net pressure and the offset distance from the pre-existing hydraulic fracture are crucial parameters that affect the FIP. This study provides an efficient approach for calculating the FIP for multi-stage HF and can aid engineers in developing tight sand reservoirs.

Thermal fracture analysis of a two-dimensional decagonal quasicrystal coating structure with interface cracks

The displacement and temperature discontinuity method is extended to the coating structure of a two-dimensional decagonal quasicrystal–crystal system, where the behavior of interfacial cracks is studied under a thermo-mechanical load. The fundamental solution of unit-point displacement and temperature discontinuity is derived on a crack surface. Based on superposition principle and finite part integrals, the boundary integral–differential equation is obtained. The algebraic equations for solving displacement and temperature discontinuities on line elements are constructed for numerical simulation. Stress intensity factors and energy release rate at crack tips are represented by the discontinuities. Influences of many factors on fracture are discussed.

Scaling Behavior of Thermally Driven Fractures in Deep Low‐Permeability Formations: A Plane Strain Model With 1‐D Heat Conduction

AbstractInjection of cold fluids through/into deep formations may cause significant cooling, thermal stress, and possible thermal fracturing. In this study, the thermal fracturing of low‐permeability formations under one‐dimensional heat conduction was investigated using a plane strain model. Dimensionless governing equations, with dimensionless fracture length , aperture , spacing , time , and effective confining stress , were derived. Solution of single thermal fracture was derived analytically, while solution of multiple fractures with constant (or dynamic) spacing were obtained using the displacement discontinuity method (and stability analysis). For single fracture, increases nonlinearly with and then transitions to scaling law , indicating that late‐time fracture length increases linearly with the square root of cooling time. For constantly spaced fractures, deviates from the single‐fracture solution at a later for a larger , showing slower propagation under inter‐fracture stress interaction. For dynamically spaced fractures, fracture arrest induced by stress interaction was determined by the stability analysis; the fully transient solution provides evolution of dimensionless fracture length, spacing, aperture, and pattern; a similar scaling law, with , obtained shows the effect of both stress interaction and fracture arrest. The solution and scaling law provide fast predictions for all reservoir and cooling conditions using (single) model parameter . Application to a geothermal site with demonstrates that thermal fractures reach 0.67, 6.25, and 78.00 m in length, 0.49, 2.30, and 13.00 m in spacing, and 0.43, 2.09, and 12.19 mm in aperture at 1, 100, and 10,000 days.

Fully 3D Simulation of Hydraulic Fracture Propagation in Naturally Fractured Reservoirs Using Displacement Discontinuity Method

Summary The shape, size, and orientation of natural fractures significantly impact the geometry of hydraulic fractures in unconventional reservoirs, such as shale gas/oil, tight gas, and enhanced geothermal system. The behaviors after the hydraulic fracture encounters natural fractures have been summarized as crossing, diverting, step over, and stopping based on a great number of numerical and experimental analysis with plain strain assumptions. However, under practical situations, the geometries of natural and hydraulic fractures are much more complex than the vertical and rectangular shape assumed by 2D models. The experimental studies on the role of height and inclination of natural fractures have revealed some other phenomena (e.g., bypassing process that the fluid-driven hydraulic fracture propagates up the back side of the natural fracture), which are unable to be captured in 2D. To better describe the intersection behaviors among different kind of fractures, a fully 3D model based on the displacement discontinuity method (DDM) is developed on top of our previous models, which only considered the propagation of hydraulic fractures. A novel crossing criterion to judge whether the hydraulic fracture will cross the cemented natural fracture in 3D is proposed. The successive node snapping scheme is adopted to construct conforming meshes for the evolving intersected curves between hydraulic and natural fracture surfaces, which only alters the location of a small fraction of nodes without changing the nodal connectivity. With this model, the evolution of fracture geometry after the hydraulic fracture intersects with the natural fractures of different toughness, size, orientation, and number is investigated. Because an extra dimension is considered, the fractures are allowed to propagate in more directions, resulting in a series of complex fracture geometries. The dynamic grid evolution method proposed in this work can promote the development of DDM in modeling fully 3D fracture networks in naturally fractured reservoirs.

Investigation of Strain Fields Generated by Hydraulic Fracturing with Analytical and Numerical Modeling of Fiber Optic Response

Summary The use of fiber optics in reservoir surveillance is bringing valuable insights into fracture geometry and fracture-hit identification, stage communication, and perforation cluster fluid distribution in many hydraulic fracturing processes. However, given the complexity associated with field data, its interpretation is a major challenge faced by engineers and geoscientists. In this work, we propose to generate distributed strain sensing (DSS)/distributed acoustic sensing (DAS) synthetic data of a crosswell fiber deployment that incorporates the physics governing hydraulic fracturing treatments. Our forward modeling can be used to add value to the interpretation task. The forward modeling is based on analytical and numerical solutions. The analytical solution is developed integrating two models: 2D fracture (e.g., Khristianovic-Geertsma-de Klerk known as KGD) and Sneddon’s induced stress. DSS is estimated using the plane strain approach that combines calculated stresses and rock properties (e.g., Young’s modulus and Poisson’s ratio). On the other hand, the numerical solution is implemented using the displacement discontinuity method (DDM), a type of boundary element method, with net pressure and/or shear stress as the boundary condition. In this case, the fiber gauge length concept is incorporated deriving displacement (i.e., DDM output) in space to obtain DSS values. In both methods, DAS is estimated by the differentiation of DSS in time. The analytical technique considers a single fracture opening and is used in a sensitivity analysis to evaluate the impact that rock/fluid parameters can promote on strain time histories. Moreover, advanced cases including multiple fractures failing in tensile or shear mode are simulated applying the numerical technique. Results indicate that our models are able to capture typical characteristics present in field data: heart-shaped pattern from a fracture approaching the fiber, stress shadow, and fracture hits. In particular, the numerical methodology captures relevant phenomenon associated with hydraulic and natural fractures interaction, which is often interpreted purely in terms of opening fractures. We can anticipate that the developed forward modeling, when embedded in a classification or regression artificial intelligence framework, will be an important tool adding substantial insights related to field fracture systems that ultimately can lead to production optimization. Also, the development of specific packages (commercial or otherwise) that explicitly model both DSS and DAS, incorporating the impact of fracture opening and slippage on strain and strain rate is still in its infancy. This paper is novel in this regard and opens up new avenues of research and applications of synthetic DAS/DSS in hydraulic fracturing processes.

Analysis of heterogeneous structures of non-separated scales using curved bridge nodes

Numerically predicting the performance of heterogeneous structures without scale separation is a challenging task owing to the critical requirements related to computational scalability and efficiency that must be satisfied. In addition, adopting a sufficiently fine mesh to consider the small-scale heterogeneities results in prohibitive computational costs, whereas neglecting them tends to drastically over-stiffen the rigidity of the structure.Thus, this study proposed an approach for constructing new material-aware shape (basis) functions per element for a coarse discretization of the structure considering each curved bridge node (CBN) that is defined along the boundaries of the elements. Rather than formulating their derivation by regarding them as a nonlinear optimization problem, the shape functions were constructed mapping the CBNs to the interior nodes and were subsequently presented in an explicit matrix form as a product of Bézier interpolation and boundary–interior transformations. The CBN shape function captures the heterogeneity of the coarse element with greater flexibility, overcomes the important and challenging issues of inter-element stiffness and displacement discontinuity across interfaces between coarse elements, and improves the analysis accuracy by several orders of magnitude. Moreover, they satisfy the basic geometric properties of shape functions, thereby avoiding non-physical analysis results. Furthermore, the performance of the proposed approach was tested and demonstrated through extensive numerical examples, including a 3D industrial example of billions of degrees of freedom, and comparisons with results obtained from classical approaches were made.